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Abstract 

This paper presents theory for the Lagrange co-rotational (CR) formulation of finite el- 
ements in the geometrically nonlinear analysis of 3D structures. In this paper strains are 
assumed to be small while the magnitude of rotations from the reference configuration is 
not restricted. A new best fit rotator and consistent spin filter are derived. Lagrange CR for- 
mulation is applied with Hybrid Trefftz Stress elements, although presented methodology 
can be applied to arbitrary problem formulation and discretization technique, f .e. finite vol- 
ume methods and lattice models, discreet element methods. Efficiency of CR formulation 
can be utilized in post-buckling stability analysis, damage and fracture mechanics, mod- 
elling of dynamic fragmentation of bodies made from quasi-brittle materials, solid fluid 
interactions and analysis of post-stressed structures, discreet body dynamics. 

Key words: Lagrange and co-rotated formulation, large rotations, foam, geometrically 
nonlinear, Trefftz elements 



1 Introduction 



Nowadays the Lagrangian kinematics description for Finite Element Method of ge- 
ometric nonlinear structures is utilized mainly in two solution techniques, i.e. Total 
Lagrangian Formulation and Updated Lagrangian Formulation. In this paper we in- 
vestigate a less common solution technique for geometrically nonlinear problems, 
i.e. the Lagrange Co-rotational (CR) formulation. 

Key features of the Lagrange CR formulation are efficiency and robustness, i.e. 
computation of tangent matrices is is low-priced, the method is easy to parallelize 
and can be used with various discretization techniques. Large class of problems 
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with finite rotations but small strains, e.g. post-buckling stability analysis, damage 
and fracture mechanics, modelling of dynamic fragmentation of bodies made from 
quasi-brittle materials, solid fluid interaction, analysis of post-stressed structures 
or discreet body dynamics can be solved efficiently with the use of Co-rotational 
formulation. 

As noted above, the Lagrange CR formulation can be applied together with dis- 
cretization techniques not suitable for Total Lagrangian Formulation or Updated 
Lagrangian Formulation. In this paper the CR formulation is applied with Hybrid 
Trefftz Stress (HTS) finite elements. Application of the CR formulation together 
with HTS elements results in a powerful synergy, which is exemplified on a case 
study problem presented in the paper. Although the presented methodology is ap- 
plied to problems of solid mechanics, the presentede method can be applied to 
image recovery, where extraction of rigid body motion from image is an issue. 

This paper follows earlier contributions BI4I5I6I1 . Full literature review can be fond 
in [O. Despite the large body of literature on co-rotational formulations, formu- 
lation of a best fit rotation functional remains an open reasearch question. In this 
paper we introduce a new general best-fit rotator functional. Additionally a new 
Spin-Filter operator, mapping increment of degrees of freedom into increment of 
rotation vector, consistent with new best-fit rotator functional is determined. 

In the first section approximation functions for the HTS elements are briefly de- 
scribed. Next, kinematics of rotating deformable motion is described and the best- 
firt rotator functional is formulated. Section|4]is supplemented with linearized equa- 
tions, used to determine the rotation operator. In Section[5]a Spin-Liver linear oper- 
ator transforming an increment of rotation angle into an increment of displacement 
degrees of freedom, and Spin-Filter are formulated. With projection matrices at 
hand, a linearization equations of HTS elements are shown in Section |4} In Section 
[7] two numerical examples are presented. The paper is concluded in Section [8} 



2 Approximation functions 



Lagrange CR formulation can be applied to a large class of discretization methods. 
In this paper we apply it to a particular instance of the Hybrid Trefftz Stress ele- 
ments. In authors view combination of the HTS element with the CR formulation 
is a good example of the robustness of the co-rotational approach: we apply the 
CR formulation to a discontinuous displacement approximation, exclusively deter- 
mined only on element boundary. Moreover, unlike for the classical displacement 
Finite Element Method, the employed displacement approximation does represent 
a partition of unity. 
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2.7 Approximation of displacements on element faces 




Fig. 1. Coordinate system on the reference element face and the element face in the refer- 
ence configuration in space 

In the spirit of hybrid FE formulations, displacement field is approximated on ele- 
ment boundary, independently on each element face. In the presented approach we 
used simple monomials given in the base coordinate system of element faces, see 
Eq. ([2]). Displacements of face / at point X/ (expressed in face coordinates) read 



u(X/) = Ur(X/)q. 
Matrix of displacement approximation functions has simple form 



(1) 



Ur(X/) = 



1 
1 
1 



Xi X2 
Xi X2 
li I2 



xl 
xl 
xl 



(2) 



where X is a coordinate vector in the reference system of element face. 



Such approximation technique allows for an easy extension of the approxima- 
tion space by arbitrary functions, which are not necessarily a partition of unity. 
Moreover, since degrees of freedom are associated with elements faces, for three- 
dimensional tetrahedrons (utilized here), an element face can have no more than six 
neighbours (adjacent through elements on both sides). As a result, this formulation 
is suitable for parallel computing where intercommunication between processors 
affects efficiency. 
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2.2 Approximation of stresses in element volume 

Detailed description of HTS elements can be found in II1I2II . In this paper we limit 
ourselves to basic equations. In this section approximation of stresses is briefly 
described. The Cauchy stress field within an element is approximated directly as: 

a = Sv(X/)v (3) 

where Sy is a matrix of field approximation functions and v is the unknown vector 
of generalised stress degrees of freedom. We note that Sv(X/) is defined on element 
face in Lagrange (co-rotated) coordinates. In eq. ([3]), the stress approximation field 
is chosen in order to automatically satisfy the equilibrium condition (??): 

L^Sv = 0. (4) 

where L is a differential operator related to the equilibrium equations. According to 
the Cauchy equilibrium condition, the corresponding stress induced traction field 
on the element boundary is 

t = na = nSyV (5) 

Furthermore, assuming small strains and elastic body, the displacement field within 
the domain is directly associated with this stress approximation expressed as: 

u = UvV. (6) 

where Uy is a displacement approximation matrix, such that Sy = C~^LUv, where 
C is a compliance matrix. 

The stress field and the associated displacement field are derived from appropriate 
Trefftz functions. For example, the stress and displacement field could be given 
as polynomial functions derived from the Airy's stress function. Alternatively, the 
Trefftz functions can be extended to functions which can take into account sin- 
gularities or are suitable for a particular boundary value problem. Moreover, the 
approximation of stress is not a priori dependent on a particular constitutive law. 



3 Kinematics of CR frame 

For completeness some basic equations and theory of rotations are presented. 
Consider a body under rotational motion 

x(X) = RX (7) 
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where R is 3x3 is the rotation operator (also called rotor), X is a selected material 
point and x is the spatial image of X. Since R is the rotation operator, distances 
between points and orientation are preserved, i.e. RR^ = I and det(R) = 1. 

In mathematical terms the set of all 3 x 3 orthogonal matrices with positive de- 
terminant forms the special orthogonal group 50(3), which is a manifold. Since 
rotations are points of the manifold 50(3), at each point R in S0{3) a tangent 
Euclidean space TrS0{3) collects rotation velocities R. 

Consider now a nine-dimensional space of all 3 x 3 matrices and a surface defined 
by the constraint /(R) = RR^ — 1 = 0. The surface normal is orthogonal to vectors 
R G T-rS0{3) or in other words 



RR^ + RR^ = 0. (8) 

Let us define co as 

CO = RR'^. (9) 

Since equation ^ is satisfied for all R G 50(3) and R G TrS0{3) we note that CO 
is anti-symmetric, co^ = -co. 

The operator co bears the name of the spatial angular velocity 

X = cox. (10) 



For CO = const the ordinary differential equation (10) has a known solution 



x(0=Exp(fco)R(0)X(0). (11) 

where time is expressed explicitly by t and Exp stands for the matrix exponential. 
Although in general CO is not constant in time, the above formula serves as a useful 
device for implementing incremental updates of R. 



3.1 Variations of rotation operator 



With the brief theory at hand, we introduce variations in order to construct tangent 
operators for the HTS formulation. From ^ we can see that R = coR and hence 
the variation of R can be expressed in the spatial form 



6R = 60R (12) 

where d^eT^S0{3). 
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3.2 Variation of displacements 



For simplicity and without loss of generality we analyse motion without the rigid 
body translations. New position of a particle is given by 



x"(X)=RX 



(13) 



where X is the position vector in the reference configuration. With that at hand we 
define displacement by 



u''(X)=x''(X)-X=(R-I)X 



(14) 



Displacements are additively decomposed into the deformational and the rotational 
part, which yields 



u^(X) = u(X) + (I - R)X = u(X) - u'*(X) 



(15) 



where and are the deformational and the rotational parts of displacements, 
respectively. 

In order to derive tangent operators, a useful variation of deformational displace- 



ments ( 15 ) is given by 



6u^(X) = 5u-5RX 



= 6u-Spin[63>]R(X-a) 
= 6u-Spin[63>]x' 
= 6u + Spin[xn6d 



(16) 



where variation of the axial rotation vector is used. We note, that 63> is a pseudo- 
vector, i.e. a 3-dimensional representation of a 3 x 3 antisymmetric matrix. The 
operator Spin[ ] takes form 



(17) 








-■^3 


X2 


Spin[x] 







-Xi 












for an arbitrary vector x. 



4 Best fit rotor 



Computation of the rotation operator for a given displacement field is at the hart of 
the co-rotational formulation. In our view, an ideal co-rotational formulation should 
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address the following objectives: 



• Versatility. It must work for static and dynamic problems, with arbitrary dis- 
cretization of displacements. 

• Rigid bodies. For rigid body motion, it should give the same results as a body- 
attached frame. This simplifies coupling of FEM-CR methods with multibody 
dynamics codes [6J. 

• Finite stretch and pure sheer If a element is subjected to a uniform finite stretch 
or pure sheer, no spurious rotation should appear. 

In order to find the rotation operator R, a new functional is proposed here. Sta- 
tionary points of this functional determine rotation, which complies with the above 
objectives. 

We note, that for a rotation free deformation, the antisymmetric part of the displace- 
ment gradient disappears. With that observation at hand, assuming that the rotation 
is constant in the element domain, the antisymmetric part of the volume averaged 
deformational displacement is used to construct the best fit rotator. 



Noting that the antisymmetric part of the displacement gradient has three non-zero 

linearly independent components, the antisymmetric part of the displacement gra- 

— * 

dient can be uniquely expressed by a pseudo- vector h. Utilizing that and applying 

— * 

integration by parts, a pseudo-vector h can be expressed by the boundary integral 



h = X u^dr = ^Spin[n]u^dr 



(19) 



where n is the spatial normal vector field on the element boundary F in the current 

or reference configuration (since the motion is roughly rigid we can assume the 

— * 

Jacobian to be nearly one). A length of the pseudo- vector h is used to formulate the 
best fit functional. Utilizing that /pSpin[N]XdF = 0, the functional takes form 



J(R) = l[h]T[h] = i 



J Spin[n]u' 
J Spin[RN](u-u'")dr 
^Spin[N]R'^xd; 



dr 



^Spin[n]u'='dr 
Spin[RN](u-u'")dr 
^Spin[N]R^xdr 



(20) 



where N is the referential normal vector field on the element boundary F in the 
current or reference configuration. We note, that the functional expressed in the 
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above form can be applied to problems where displacements are only known on 
elements boundaries (and can be discontinuous). 



For given displacements u, a new position vector x in the current configuration is 
easy to determine. For a fixed x, the condition of stationarity of J(R) represents the 
set of nonlinear equations for finding R 



dR 



^Spin[N]R^ 
J Spin[N]R^ 



xdr 



xdr 



^Spin[N]5R^xdr 
J Spin[N]R'^Spin[x]dr 



5^ = 0. 



where ( |T2| ) has been utihzed. We then requre that 

T r 



H 



^ SpinfNlR'^xdr ^ Spin [NjR^^ Spin [x]dr 







and solve the above sytem by means of the Newton-Raphson iterations 



an 



-A<I> 



H 



3$ 

Rk+i = Exp(A$)Rk. 



In order to approximate (23 ) we linearise Exp(A<I>) 



^ /.^N , sin||A<I>|l ^ 1 -cos II A<I> 11 ^2 
Exp (A<I>) = I + .. " ^ .. " A<I> + J'... " A^^ 



A4) 



A4)||2 



I + A$ 



(21) 



(22) 



(23) 
(24) 



(25) 



and noting that Rk+i f» Rk(I + A4>) we represent (23 ) by up to first order terms of 
H(Rk(I + A<I>)) =0. That is 



^Spin[N]RjSpin[x]dr 
J Spin[N]Rjxdr 



Spin[N]RjSpin[x]dr 



A4) 



Spin[N]Rjxdr 



J Spin[N]RjSpin [Spiny[x]]dr 
Spin[N]RjSpin[x]dr 



(26) 



where e,- is the Cartesian base vector, i,j G {1,2,3} and Spin^ [x] represents the yth 
column of Spin[x]. 
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In this section a new best-fit functional, based on boundary integrals is presented. 
Stationary point of this functional lead to a system of non-linear algebraic equa- 
tions, form which current rotation operator can be calculated by means of Newton 
iterations. It should be noted, that the above procedure allows to determine the rota- 
tion quite independently of the motion history. However, in order to avoid troubles 
with multiple minima and correctly trace the history of rotation it is appropriate to 
initialize the iterative procedure with the most recent value of rotation. 



5 Projection matrices 



In order to construct tangent operators for the consistent Newton scheme projection 
matrices filtering a variation of the deformational displacements from the variation 
of the total displacements have to be derived. 

5q^ = P5q = (I-SG)6q (27) 



Projection matrix P, Spin-Liver S and Spin-Filter G depend on the approximation 
method. In this paper we construct these matrices for the particular approximation 
of HTS finite elements. However, the methodology is general and can be applied to 
a large class approximation methods. 



5.0.7 Spin Liver 

A Spin-Liver matrix represents a linear operator transforming a variation of the 
rotation axis into a variation of linear displacements. Spin-Liver S depends only on 
constant and linear approximation functions. Higher order terms in displacements 
approximation basis are neglected in the formulation of Spin-Liver. It can be noted, 
that for any displacement function given by higher-order terms the motion is not 
stress free, hence higher-order therms cannot be a part of the rigid body motion. 

Exploiting features of the HTS displacement approximation, see Eq. ([2]), the Spin- 
Liver S can be defined for each face independently, which yields 

6u^(X) = Ur(X/)S/(X^sc)5^ (28) 
where X^^^ = X(X = 0) is a face geometrical center. The Spin-Liver matrix itself, 
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Fig. 2. For a given current face frame and an arbitrary variation of spin 8fl, the variation 
of rotational linear displacements are presented in figure on the left-hand side. Linear rota- 
tional displacements are computed using formula 8u^ = S55>. The green triangular face rep- 
resents the reference configuration, the red triangular face depicts the current configuration. 
The Spin-Liver is constructed so that the resulting variations of the linear displacements are 
orthogonal to the variation of the rotation axis and tangent to the sphere, see figure on the 
right-hand side. Position of points on the sphere depends on the current rotation axis and 
rotation angle. 

for our particular implementation of the HTS element has simple form 






4 


Xj 











-A 








3x5 /3Xi 


-3x^/3Xi 







3x;/3Xi 




-3x^/3Xi 








3x5/3X2 


-3x^3X2 







3x^/3X2 


3x^3X2 


-3x^/3X2 















vv^here = RX^^^ and 3x^/3X|xfgc = R3X/3X|xfgc. Definition of the Spin-Liver for 
other types of elements can be found in 
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5.0.2 Spin Filter 




Fig. 3. Current configuration is given by the green tetrahedron. Variation of current con- 
figuration is depicted by the blue tetrahedron. Variation of the rotation axis is visuaHzed by 
the black arrow. Variation of the deformation is depicted by the red tetrahedral. 



Matrix G represents a linear operator which transforms a variation of linear dis- 
placements into a variation of rotation axis. In other words Spin-Filter accounts for 
a change of magnitude and direction of the rotation axis due to the change of linear 
displacements. 

In order to get quadratic convergence, Spin-Filter matrix has to be consistent with 
functional /, defined in Section |4} For a given current configuration, i.e. rotation R, 
and linearized displacements, functional / takes form 



/(63>) 

_ 1 
~ 2 



-| r r 

J Spin[N]R'^(6u + Spin[x'']6d)dr J Spin[N]R'^(6u + Spin[x'']63>)dr 

(30) 
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Stationary point of this functional yields 



36^ 



Spin[N]R'^Spin[x'"]dr 



Spin[N]R'^(6u + Spin[x'"]63>)dr 



(31) 



Equation of the stationary point is given by the system of three linear equations 

T 



6<I> = 



J Spin[N]R'^Spin[x'"]dr J Spin[N]R'^Spin[x'']dr 

r 1 T r 

^ Spin[N]R'^Spin[x'"]dr ^ Spin[N]R'^6udr 

Substituting 5u = Ur5q yields 



(32) 



J SpinfNjR'^Spinfx^dr J Spin[N]R'^Spin[x'"]dr 

^ Spin[N]R'^Spin[x'"]dr ^ Spin[N]R^Urdr 



(33) 



5q 



As a result, a matrix representing the Spin-Filter 

= G6q, 

is expressed in the form 

T 



(34) 



G = 



^ Spin[N]R^Spin[x'"]dr ^ Spin[N]R'^Spin[x'"]dr 

J Spin[N]R^ 



(35) 



Spin[N]R'^Spin[x'"]dr 



Urdr 



The novelty of this approach in comparison to those found in other papers related 
to the CR formulation, is that the Spin-Filter is derived form best fit functional and 
can be therefore applied to a large class of approximation methods. 



6 HTS formulation 



Although the co-rotational formulation can be derived independently form a prob- 
lem formulation by consistent differentiation of projection matrices [6J, in this pa- 
per we will formulate tangent operators and residuals directly from weak forms. In 
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reference current corotated 



deformed element 




Fig. 4. Refrence and coorotated frames 



authors opining this allows for better understanding of the problem, without loss of 
generality. 

The co-rotational Hybrid Trefftz Stress element is built using Lagrange coordinates, 
i.e. integrals are over reference volumes and faces. Stresses are expressed in the co- 
rotated (current) configuration and displacements are expressed in the reference 
coordinate system. At each iteration co-rotated (Cauchy) stress is pulled-back to 
the reference coordinate system (or conversly, the displacement is pushed-forward 
to the co-rotated reference system). 



6,1 Integral form 



Static boundary conditions are enforced in strong integral form, i.e. we weight static 
boundary condition by function wi and integrate over element boundary. 



Function wi is given in the reference coordinate system and stresses are pulled- 
back to the reference coordinate system. This results in geometrical nonlinearities, 
i.e. the integral becomes a function of the unknown rotor R. 




(36) 
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6,2 Weak form 



The compatibly conditions within a body domain are imposed in a weak sense, 
where by the differentiation by parts Dirichlet boundary conditions are included 

tr[w28]d^2- / (W2n)^(ur-c-u'')dr- / (w2n)^udr = (37) 
[ tr[w28]d^2- / (W2n)'^ufdr- / (w2n)'^udr = (38) 

Ja JTu JTa 

Re= / tr[Rw2R^8]d^2- / (Rw2ii)^ufdr- / (Rw2n)^udr = (39) 
Ja Jr^ Jvcy 

Re= f tr[Rw2R^ReR^]d^2- [ (Rw2n)^ufdr- / (Rw2n)^udr = (40) 

f tr[w28]d^2- [ (w2n)^R^ufdr- / (w2n)^R^udr = (41) 

A weighting function W2 is defined in the co-rotated coordinate system and is trans- 
formed to the reference coordinate system. We note, that the term under the volume 
integral is a scalar, so it not depend on the coordinate system. As a result it is con- 
stant for all steps analysis and is computed only once. 



6. 3 Linearization 



The system of nonlinear algebraic equations is linearized and solved using Newton 
method. The rotation operator R is a function of current displacements and it is 
computed at each iteration using the best-fit functional ( [20| ). At a current cofigura- 
tiom, linearized Eq. ( |36| ) takes from 

Ra + ^ wi'^Reandr + wi'^eRdndr = (42) 

Ra + y^ wi'^R6andr + j wi^6<I>Randr = (43) 
Ra + wi^RSandr - j wi^Spin[Ran]dr63) = (44) 

In a similar manner to Eq. ([36]), residual Eq. ( |4T] ) is linearized and yields 

Re + y^ tr[w258]d^2 - y^(w2n)^R^5uf dr - y^(w2n)^6R^uf dr = (45) 

Re+ / tr[w268]d^2- / (w2ii)^R^6ufdr+ / (w2n)^R^5cI>uf dr = (46) 
Ja Jr Jr 
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Re + y^ tT[w2bE\d^l- J (w2n)V6ufdr-y^ (w2n)'^R'^Spin[uf ]dr6d = 

(47) 



A key term in above equations is the change of the angular vector 50. We note 
50 depends on the current configuration and is a hnear function of the variation of 
displacement degrees of freedom, see Eq. ([34]) and Eq. ( [35] ). 

6.4 Matrixform 



The linearized equations are reformulated, with stresses and displacements sub- 
stituted by suitable products approximation functions and increments of vectors 
storing stress and displacement degrees of freedom. Spin-Filter G and Spin-Liver S 
are utilized in order to express linearized equations exclusively in terms unknown 
vectors, which yields 

Ra + ^Uj^nSydrSv - ^Uj;Spin[Ran]drG6q = 0, (48) 

Re + y^ iiSv^Uvdr6v-y nSv^R^UrdrSGSq-y iiSv^R^Spin[uf]drG6q = 0. 

(49) 

The above linearized equations can are expressed conveniently in the matrix form 



F -AT(I-SG)-QG 
-A UG 





6v 




-Re 











where 



r = ^nSv^Uvdr, 
A = ^nSv^R'^Urdr, 
Q = ^nSv^R^Spin[uf]dr, 
U = ^UfSpin[Rdn]dr. 



(50) 

(51) 
(52) 
(53) 
(54) 



6,5 Partial Solution 



The size of the problem to be solved can be reduced by static condensation, taking 
advantage of the element level approximations of stress fields. On element level, 
solution for stress degrees of freedom is expressed by 

5v = -F-^Re + F-^ { aT(I - SG) + QG} 6q. (55) 
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Substituting 5v in Eq. ( [48] ) yields 

{-AF^ [aT(I-SG)+QG] -UG}6q = -Ra-AF^Re. 



(56) 



This leads to a modified system of equations 



F6q=-Ra. 



(57) 



In the above, the stress degrees of freedom are eliminated conveniently from the 
global system of equations, leaving only the displacement degrees of freedom to be 
determined. This does not affect the final result in any way and simply represents 
a solution technique that reduces the number of equations that have to be solved 
simultaneously. Following the solution of the displacement degrees of freedom, 
the stress degrees of freedom can be recovered on an element by element basis, a 
process that can be easily parallelized. 

It is noted, that independent displacement approximation on elements faces, as pre- 
sented here, lead to larger number of unknowns in comparison to classical FEM. 
Nevertheless, the matrix sparsity, which has a significant impact on the efficiency 
of direct and iterative solvers, is minimised in the presented approach due to the 
association of the displacement degrees of freedom with faces rather than vertices 
and due to the fact that any given face can only have two neighbouring elements in 
either 2D or 3D. For the same reason the communication between parallel proces- 
sors is significantly reduced. As shown in [?], this approach results in very good 
scalability. 



7 Numerical examples 

To illustrate features of the CR-HTS approach, two numerical examples are pre- 
sented. The first one is an academic problem investigating model's accuracy. The 
second example further illustrates model's accuracy for large scale 3D representa- 
tive model of complex micro structure. 

Examples are analysed with the use Portable Extensible Toolkit for Scientific Com- 
putation (PETSc) [7J and A Mesh-Oriented datABase (MOAB) [SJ. Finite elements 
and CR formulation were implemented in C and C++. 

7.7 Pure bending in 3D 

In the first numerical example the numerical results are compared with an analytical 
solution. A beam (0.2x0.2x5.0) with square cross-section is stressed under pure 
bending. 
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In order to get analytical solution, a special case is analysed, i.e. E = I and v = 0. 
A plane curvature is given explicitly as 



K = ' (58) 

(l+w'2)3/2 

and the bending moment according to Bernoulli theory is given by 

M = ElK. (59) 

Strain energy is calculated from 



1 



ElK^d/ (60) 



1/2 



Fig. 5. Mesh 




Fig. 6. Deformation for angles of rotation k and 2k at the end of the bar. The deformation is 
plotted using displacements Trefftz approximation functions (Uv) supplemented with rigid 
body motion. 

With the analytical solution at hand, a 3D numerical model is analysed. A dis- 
cretized beam consisting 1 160 tetrahedrals as shown in Fig.[5| A HTS element with 
quadratic displacements discretization and 3rd order polynomial approximation of 
stress field within element is used. At both ends the kinematic boundary conditions 
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Fig. 7. Comparison of the analytical and numerical solutions, 
are applied in order to force pure bending 



0, V = }^cos^Lk — _y, 



K 



. 1 

sm -Lk — z 



(61) 



where the length of the beam L = 5. Despite the coarse mesh, numerical solution 
is in excellent agreement with the analytical one, cf. Fig. |6| For K = 27r/L beam 
forms a circle, which is in agreement with the analytical solution. At each time step 
quadratic convergence is observed. 

This test has been reproduced for several bar placements in 3D, in order to verify 
the implementation and the formulation of best-fit rotator. We note that the Cauchy 
stress field distribution is in an exact much which the analytical solutions for all 
cases. 

At each increment an updated rotor R is computed for each element independently. 
At all times the procedure presented in section |4] required few iterations to achieve 
the machine precision. Moreover, profiling the runtimes showed that the time re- 
lated to the rotor computations was negligible. 



7.2 Analysis of foam 



In order to present one of the many potential applications of the Lagrange CR 
formulation we investigate deformation of a hypothetical foam. The mesh was ob- 
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Fig. 8. Foam surface mesh and geometry 



tained from a high-resolution 3D scan, included as a case study in the Simpleware 
software The mesh comprised 95956 tetrahedrons and 40725 nodes. The con- 
stitutive parameters of the model were chosen arbitrarily, i.e. E = le4 and v = 0.25. 
At the bottom plate a fixed boundary conditions were applied. At the top plate trac- 
tions in z direction were applied. On the front, back, left and right side the displace- 
ment in the normal direction to the side wall was blocked. 

Displacement and stresses filed were approximated using Hybrid- Trefftz elements. 
Two cases were considered with linear and quadratic approximation of displace- 
ments. For the first case, i.e. linear functions approximating displacements, poly- 
nomial Trefftz functions of order two were used to approximate stresses. For the 
second case, polynomial Trefftz functions of order three were used to approximate 
the stress field. This led to the model with ^ 5 million DoFs and ^10 million 
DoFs, for the case one and case two, respectively. We note, that we not observe 
shear locking phenomenon for the HTS elements LlOl . despite tetrahedrons being 
used to approximate the geometry. 

Large rotations led to local buckling in the micro-structure which had an impact 
on the convergence of the Newton method. In order to improve efficiency, a set of 
nonlinear equations was supplemented with the arc-length control. We introduced 
additional unknown X, such that 



R<7 = Ra — ^T;^5 



(62) 



and additional control equation 



R;, = GAq^GAq + ^^^i^ - s^. 



(63) 
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GAq is an increment of the rotation angle AO, which yields 



(64) 



where elements of matrix G are determined at the beginning of a load increment. 
^ and X are scaling parameter and load factor respectively, s is radius of a hyper- 
sphere. We note that the change of the rotation angle is a source of nonlinearities 
and a direct control of the rotation angle positively affects the stability of the New- 
ton method. Since controlling geometrical instabilities in a structure like foam is an 
open problem, the above solution has been adopted only as a demonstration, while 
further improvements could well be introduced. 
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Case 1 , linear approximation of displacements 

Case 2, quadratic appro:Jcimation of displacements 

Linear response 

O snapshot 




Fig. 9. Equilibrium path 

Analysis of foam deformation under compression and tension showed expected 
foam response, i.e. overall response softening and stiffening for compression and 
tension, respectively. Deformation, rotation angle and stress Cxx are shown in Fig. [10 
and Fig. [TTJ We can note that for compression we observed interpenetration and 
the model should be extended to take into account self-contact. For compression, a 
transition to bending dominated behavior was observed. For tension, a straighten- 
ing of foam's branches was observed. In both cases local buckling was a source of 
numerical difficulties. 



Comparison between linear and quadratic approximation of displacements, i.e. case 
1 and case 2, showed a close much between the numerical responses. Good quasi- 
quadratic convergence for case 1 and case 2 was achieved until the final point of 
analysis. However poor mesh quality and local instabilities were the source of prob- 
lems with the convergence for extreme (non-physical) loads. Those problems need 
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further investigation and in authors opinion can be overcome by a better local arc- 
length control and line searches. 




Fig. 10. Rotations, see Fig. |9] where snapshots are taken. 




Fig. 11. Stress o^x, see Fig. [9] where snapshots are taken. 

In figures showing deformation, displacements are plotted using Trefftz displace- 
ment function u, see Eq. |6} supplemented with rigid body rotation u^. Trefftz dis- 
placement field is adjoint to stresses calculated to fulfill equilibrium, hence dis- 
placement vector u in Eq. [6] do not carry any information about translations and 
rotations. By visual inspection of the figures in this example and in Fig. |6] we can 
verify an excellent consistency of the calculated rotor with the overall deformation. 



8 Conclusions 

This paper presents a geometrically nonlinear Lagrange co-rotational formulation 
applied to Hybrid Stress Trefftz Elements. All terms in the internal force and tan- 
gent matrices are accounted for. 

A new best-fit rotator for continuum elements is developed. The best-fit rotator 
functional is constructed by utilizing a boundary integral over a deformable el- 
ement/body. Stationary point of the functional lead to the set of three nonlinear 
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equations, solved by means of the Newton method. Formulation of best-fit rotator, 
load vector and tangent matrices was verified by two examples: pure bending test 
and large deformation of foam. 

Findings from foam analysis can be used to model auxetic and non-auxetic poly- 
meric foams [fTT1| . remodelling and growth of tubercular bone [fT2l| and fracture of 
ceramic/quasi-brittle foams. 

In authors opinion exciting research opportunities arise form this contribution. The 
CR-HTS method can be applied to many problems involving quasi-brittle materi- 
als and large deformations. We note, that for fracturing materials very often large 
rotations are observed for fractured parts. Moreover the current approach can be 
applied to many dynamic problems involving large rotations and small deforma- 
tions. For models involving in addition large strains, it is possible to use the CR 
formulation together with Total Lagrangian or Updated Lagrangian formulations, 
where the CR frame could be attached to rotating parts of the model or individual 
elements. 

Another interesting research opportunity is application of the Lagrange CR formu- 
lation to very large problems, solved with Krylov iterative solver without explicit 
formation of the tangent matrix. Matrix free methodology is attractive for memory 
demanding problems. We note, that in case of the presented methodology, all inte- 
grals can be computed only at the beginning of an analysis, and the tangent matrix 
can be calculated with the use of projection matrices only. As a result, the matrix 
free methodology can be used efficiently with nonlinear problems and higher-order 
elements without the need for integration at each iterative solver step. 
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